PRODUCTS OF RANDOM MATRICES AND GENERALISED 
QUANTUM POINT SCATTERERS 
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Abstract. To every product of 2 X 2 matrices, there corresponds a one- 
dimensional Schrodinger equation whose potential consists of generalised point 
scatterers. Products of random matrices are obtained by making these interac- 
tions and their positions random. We exhibit a simple one-dimensional quan- 
tum model corresponding to the most general product of matrices in SL (2, R). 
We use this correspondence to find new examples of products of random ma- 
trices for which the invariant measure can be expressed in simple analytical 
terms. 



1. Introduction 

Products of random 2x2 matrices arise in many physical contexts: in the 
study of random spin chains, or when calculating the distribution of the natural 
frequencies of a classical random spring chain, or more generally when considering 
the propagation of a wave in a one-dimensional disordered medium [5J [3TJ 152"] . 
It is often the case that, in the presence of disorder (i.e. randomness), the waves 
become sharply localised in space. This physical phenomenon is known as Anderson 
localisation; one of its mathematical manifestations is the exponential growth of the 
product of random matrices. 

The rate of growth is called the Lyapunov exponent; it often has a physical 
interpretation in terms of the exponential decay of the transmission probability as 
the size of the disordered region grows. One method for calculating the Lyapunov 
exponent is based on a general theory developed by Furstenberg and others [6|l9ll20j. 
This method requires the explicit knowledge of a certain measure on the projective 
space, invariant under the action of the matrices in the product. Examples of 
products of random matrices for which this invariant measure can be obtained in 
analytical form are, however, very few; see for instance [6j [10l [27l [33] and the 
references therein. 

The calculation of the Lyapunov exponent need not always make use of this 
invariant measure. There are alternative approaches; see for instance [32 1 I35 [ 155]. 
Nevertheless, the problem of determining the invariant measure is interesting in 
itself, and the present paper will focus on the presentation of new explicit examples 
from corresponding examples of exactly solvable models of one-dimensional disor- 
dered systems with point scatterers [5] [31 [35] ■ In our context, the phrase "exactly 
solvable" means that the calculation of the Lyapunov exponent associated with the 
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disordered system is reduced to a problem of quadrature. Some of the models were 
solved by Nieuwenhuizen |35| (without the use of the invariant measure); some of 
them are, apparently, new. Although the work reported here is, for the most part, 
mathematically driven, these new models are of independent physical interest. To 
the best of our knowledge, all the explicit formulae for the invariant measures con- 
stitute new results. In the remainder of this introductory section, we review some 
relevant concepts and some known facts, summarise our main results, and give a 
sketch of our approach. 

1.1. Products of random matrices. Let 

Ai, A 2 , A 3 , ■■■ 

denote independent, identically-distributed 2x2 matrices with unit determinant, 
let fi be their common distribution, and consider the product 

(1.1) n n := A n A n - X ■ ■ ■ A x . 
The number 

(1.2) 7 , := ,i m 'yUM 

n— >-oo Tl 

where | • | denotes the norm on matrices induced by the euclidean norm on vectors, 
also denoted | • |, is called the Lyapunov exponent of the product. 

The product grows if the angle between the columns decreases or, equivalently, 
if the columns tend to align along some common direction. In precise mathematical 
terms, a direction in M. d is a straight line through the origin, and the set of all direc- 
tions is, by definition, the projective space P \Mr\ . The case d = 2 is particularly 
simple: any direction 

\ ( X ) : A G 



is characterised by the reciprocal, say 

z=-eR := IU {oc}, 
V 

of its slope. So we can identify P (R 2 ) with K. The calculation of the Lyapunov 
exponent is often based on the formula |H1 [20] : 

(1.3) <y»=[v(dz)[ ^(dA)ln^ 1 

JR JSL(2,R) / z 

u 

In this expression, /j, is the known common distribution of the matrices A n in 
the product, whereas v is the — a priori unknown — probability measure on the 
projective line which is invariant under the action of matrices drawn from /i. Here, 
invariance means that if 

a b 



A 



c d 



is a /i-distributed random matrix and z is a ^-distributed random direction, then 
the direction 

. az + b 
(1.4) s/(z) 



cz 
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— of the vector obtained after A has multiplied a vector of direction z — is also 
^-distributed. In the particular case where v has a density, i.e. 

it may be shown that 

(1.5) f(z) = I M (cL4) (/ o ^) (z) ^(z) . 

JSL(2,R) az 

However, there is no systematic method for solving this integral equation. 

1.2. The particular products of random matrices considered. To describe 
them, let us first remark that every A <G SL (2, R) has a unique Iwasawa decompo- 
sition 



(1.6) 



cos9 -sm9\ fe w \ A u 
sin (9 cos (9 ) \ e~ w j \ 1 



for some 0, it, io e R. This follows easily by applying the familiar Gram-Schmidt 
algorithm to the columns of A. The three parameters in this decomposition have 
simple geometrical meanings: —9 is the angle that the first column of A makes with 
the horizontal axis, e w is its magnitude, and u is related to the angle between the 
columns; in particular, u = if and only if the columns are orthogonal. 

Now, suppose that these three parameters are independent random variables. 
We use the notation 

v ~ Exp(r) 

to indicate that v is a random variable with an exponential distribution of parameter 
r, i.e. its density is given by 

re"™ l( ,oo)(w) , 

where for every set A C R, 



U(x) = 



for x G A 
otherwise 



Also, 5 X will denote the discrete probability distribution on R with all the mass 
at x. We shall provide an explicit formula for the ^-invariant measure of the 
product n„ when the matrices are independent draws from the distribution \i of A 
corresponding to either 

(1.7) 9 <~ Exp(p) , ±u ~ Exp(q) , w ~ Sq . 



(1.8) 9 <~ Exp(p) , u <~ <5q , ±w ~ Exp(q) . 

We shall also look at other closely related products: for instance, products in- 
volving matrices of the form 

_ /cosh6» sinh(9\ fc w \ A u 
\sm\i9 cosh9j \0 c" 1 " ) \Q 1 

and we shall exhibit invariant measures for such cases too. 
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1.3. The Schrodinger equation with a random potential. Our approach to 
computing the invariant measure will not make explicit use of the integral equation 



(1.5 1. Instead, we shall exploit the fact that these products arise when solving the 
Schrodinger equation (in units such that h = 2m = 1) 

(1.9) -ip"(x) + V(x)ip(x) = Eip(x) 

for a given energy E and a potential function V that vanishes everywhere except 
on a countable set of points {xj}. Physically speaking, one can think of ip as the 
wave function of a quantum particle in a crystal with impurities; the effect of the 
impurity located at Xj is modelled by the boundary condition 

where Bj £ SL(2,R). The potential V is therefore a sum of simpler potentials, one 
for each pair (xj ,Bj), known variously as point scatterers, generalised contact scat- 



terers or pointlike scatterers [5J [31 EJ E3 H3 EE] • The case (1.7) corresponds to the 
disordered version of the familiar Kronig-Penney model [30] considered by Frisch 
& Lloyd [T5J and Kotani [35] . The case (|1.8[) corresponds to a "supersymmetric 



version" of the same model, in which the Schrodinger operator factorises as 



d 2 „, . d 2 



(L10) -d^ + V{x) = ~&x^ + W{x) W ' {x) 



■l x+ W(x) 



+ W(x) 

ax 



and the superpotential W is of the Kronig-Penney type. Such a supersymmetric 
Hamiltonian is related to the square of a Dirac operator with a random mass W — a 
model that is of independent interest in many contexts of condensed matter physics 

[TIES [Ml EI]. 

The strategy for calculating v is based on the observation that it is also the 
stationary distribution of a certain Markov process {z(x)}, where 

is the Riccati variable associated with the Schrodinger equation. In the particular 
case where 

Xj+i - Xj ~ Exp(p) 

and the Bj are independent and identically distributed random variables in SL(2, R), 
one can, following Frisch & Lloyd [TI5] , show that the density of the stationary dis- 



tribution satisfies a certain integro-diffcrential equation. The cases (1.7| and (1.8) 
share a special feature: the distribution of the Bj is such that the integro-differential 
equation may be reduced to a differential equation. Furthermore, this differential 
equation is simple enough to admit an exact solution in terms of elementary func- 
tions. 

The idea of using the Riccati variable to study disordered systems goes back to 
Frisch & Lloyd [15]. The well-known "phase formalism" introduced in [H [3T] is 
another version of the same idea. The trick that allows one to express the equation 
for the stationary distribution of the Riccati variable in a purely differential form is 



borrowed from Nicuwcnhuizcn's work |35j on the particular case ( 1.7), in which the 



Dyson-Schmidt method is used to compute the Lyapunov exponent directly from a 
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so-called characteristic function. The same trick has been used by others in various 
contexts [HI [5TJ [33J . The key fact is that the density of the exponential distribution 
satisfies a linear differential equation with constant coefficients. Our results on 
products of matrices therefore admit a number of extensions; for instance when ±v 
(or ±w) has, say, a gamma or a Laplace (i.e. piecewise exponential) distribution. 
One difficulty that arises with these distributions is that the differential equation 
for the invariant density is then of second or higher order. This makes it harder 
to identify the relevant solution; furthermore, this solution is seldom expressible 
in terms of elementary functions. Without aiming at an exhaustive treatment, we 
shall have occasion to illustrate some of these technical difficulties. 

1 .4. Outline of the paper. The remainder of the paper is as follows: in Ej2j we 
review the concept of point scatterer. The Frisch-Lloyd equation for the stationary 
density of the Riccati variable is derived in Sj3j In fj4| we study particular choices of 
random point scatterers for which the Frisch-Lloyd equation can be reduced to a 
purely differential form. We can solve this equation in some cases and these results 
are then translated in terms of invariant measures for products of random matrices. 
Some possible extensions of our results are discussed in ^5] We end the paper with 
a few concluding remarks in Sj6[ 

2. Point scatterers 

Let n€i and let 5 denote the Dirac delta. The Schrodinger equation with the 
potential 

V(x) = uS(x) 
can be expressed in the equivalent form 

(2.1) -tp" = Eip, x^O, 
and 

(2.2) ^(0+)=^(0-), ^'(0+)=^'(Q-)+uip(0-). 

This familiar "delta scatterer" is a convenient idealisation for a short-range, highly 
localised potential. 

A (mathematically) natural generalisation of this scatterer is obtained when the 



boundary condition (2.2) is replaced by 



^(0+) J V no- 
where B is some 2x2 matrix. We shall refer to B as the "boundary matrix". 
In order to ascertain what boundary matrices yield a Schrodinger operator with a 
self-adjoint extension, we start with the observation that the probability current 
associated with the wavefunction is proportional to 

'0 -1\ (^'{xf 



i o ; U( 



x 



where the bar denotes complex conjugation. The requirement that the probabil- 
ity current should be the same on both sides of the scatterer translates into the 
following condition on B [TT] : 



1 \i 
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where the dagger denotes hermitian transposition. Equivalcntly, 

b\\b 22 - b 2 ibi2 = 1 and Im(&n6 2 i) = Im(&22&i2) = . 
It is easily seen that this forces [38] 

e~ ix B e SL(2,K) 

for some real number \- As discussed in Appendix[XJ for the purposes of this paper 
there is no loss of generality in setting x = and restricting our attention to the 
case of real boundary matrices. 
We write 

(2.4) V{x)=a B (x) 

for the potential with these properties, and call it a point scatterer (at the origin) or, 
as it is also known, a generalised contact scatterer or pointlike scatterer [21I51 H71I5B] . 
We remark that the Riccati variable z — ip' /ip of the Schrodinger equation with 
this potential satisfies 

(2.5) z' = -(E + z 2 ) 7 x^O, 
and 

(2.6) 2(0+) =m(z(o-)) 

where 38 is the linear fractional transformation associated with the matrix B: 

(2.7) 38(z) = . 

' W b 2l Z + b 22 

The fact that B e SL(2,R) ensures that 38 is invertible. 

In order to gain some insight into the possible physical significance of th e bo und- 
ary matrix B, we set E — k 2 , k > 0, and look for solutions of Equations (2.1) and 
([273]) of the form 



(2.8) iP(x) 



for x < 

fl in g-ifcs £ Qr x > q 



By definition, the scattering matrix S relates the incoming amplitudes to the out- 
going amplitudes via 



a_ 



Oj. / \ a" 



(2.9) ( :- ut ) = S 

Hence 
S = 



b 2 \k 2 + ifc(6 n + 622) - &12 

^2ifc 2 -ifc(6 2 2-&ii) + &i2 2ife 

2ifc(bn6 2 2 - bi 2 b 2 i) b 2X k 2 + ik(b 22 - b\x) + b 12 

The relationship between boundary and scattering matrices is discussed at greater 
length in Appendix [A] 



Example 1. For the delta scatterer defined by (|2.2|, 
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The wave function is continuous at the origin, but its derivative experiences a jump 
proportional to the value of the wave function there. We have 



S = 



1 



2ifc 



u 2ik 
2ifc u 



and £§(z) = z + u . 



The fact that 



lim S = —I , 

u— > ±00 



where / is the identity matrix, indicates that the limiting case of an infinitely large 
"impurity strength" u corresponds to imposing a Dirichlct boundary condition at 
the scatterer's position. 



B = 



Example 2. The "delta-prime" scatterer (see for instance [21 [3B]) is defined by 

where v G R. Now it is the derivative of the wave function that is continuous at 
the origin, and the wave function that jumps: 

ip(0+) -V'(O-) =vt/]'(Q). 

We emphasise that, in spite of its (widely used) name, the delta-prime scatterer 
does not correspond to using the distributional derivative 5' as a potential [2] • 



We have 



The fact that 



S 



1 



vk 2i 



2i + vk V 2i vk 



and 38{z) 



lim S 

v— >±oo 



-I 



indicates that a Neumann boundary condition is obtained in the limit of infinite 
strength v. 

The question of the possible physical significance of the delta-prime scatterer 
was considered by Cheon & Shigehara |12) . who showed that it can in principle be 
"realised" by taking an appropriate limit of three neighbouring delta scatterers. 



Example 3. Let w G R and 
(2.10) 



B = 







In this case, the scatterer produces a discontinuity in both the wave function and 
its derivative. As pointed out in [13], the Schrodinger equation (2.1) can be recast 
as the first-order system 

-i/j' -w%l) = k4> 

4/ - W(j) = kiP 

with 

W(x) — w5(x) . 

The meaning of these equations becomes clear if we introduce an integrating factor: 
d 



d.r 
_d_ 
d.c 



exp 
exp 



W(y)dy\ ^ 

/ 

W(y)dy] 4, 



k exp 
k exp 



W(y)dyj 4 
W(y) dy] V • 
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We call this scatterer the super ■symmetric scatterer. We have 
/ tanh w sech w 



S = 



\ sech w — tanh w 



and 3§{z) 



e 2w z. 



Hence the scattering is independent of the wave number k — a property consis- 
tent with the observation, made in Albeverio et al[2 , that diagonal matrices (are 
the only matrices in SL(2, M) that) yield boundary conditions invariant under the 
scaling 

ip(x) i-» \f\ip(\x) , A>0. 

However, in contrast with the previous examples, the scattering is asymmetric, i.e. 
not invariant under the transformation x i— > —x. The limit of infinite strength w 
has a clear interpretation: it corresponds to a Neumann boundary condition on the 
left of the barrier, and to a Dirichlet condition on the right. 



Example 4. Let 



B 







This point scatterer can be thought of as two neighbouring scatterers — a super- 
symmetric scatterer of strength w on the right, and a delta scatterer of strength u 
on the left — in the limit as the distance e separating them tends to 0; see Figure 
[T] For want of a better name, we shall refer to it as the double impurity. 



w 



Figure 1. The double impurity: the empty dot corresponds to 
the location of a supersymmetric scatterer of strength w while the 
black dot corresponds to the location of a delta scatterer of 
strength u. 



We have 

S = 



1 

2ik cosh w — ue u 



and 



2ifc sinh w + ue u 
2ik 



2ifc 

-2ifc sinh w 



,%(z)=e 2w (z + u). 

This particular scatterer is interesting for the following reason: the Iwasawa 



decomposition (1.6) implies that any point scatterer for a real boundary matrix 
can be thought of as a double impurity "up to a rotation". For example, the 
boundary matrix for the delta-prime scatterer may be decomposed as 



cosf 
sin t 



— sm( 
cos 9 







PRODUCTS OF RANDOM MATRICES AND QUANTUM POINT SCATTERERS 



9 



with 



arctanu, w = — ln(l + v 2 ) and u 



We shall return to this point in the next section. 

3. A GENERALISED KRONIG-PENNEY MODEL WITH DISORDER 

In this section, we elaborate the correspondence between disordered systems 
with point scatterers and products of random matrices. Then, for a particular 
type of disorder, we show how, following Frisch & Lloyd 19], one can derive a 
useful equation for the stationary density of the Riccati variable associated with 
the system. 

3.1. The generalised Kronig Penney model. Given a sequence 

{Bj} C SL(2,R) 



and an increasing sequence {xj } of non-negative numbers, we call the equation (1.9) 
with the potential 

oo 

(3.1) \ I ■'• : x i) 

a generalised Kronig-Penney model. The notation o~b(x) was defined in Equa- 



tion (2.4 1 



Let us consider first the case where the energy is positive, i.e. E = k 2 , k > 0. In 
principle, one could dispense with the parameter k and set it to unity by rescaling x 
but, as we shall see later in j ]3.5[ there is some advantage in making the dependence 
on the energy explicit. For xj < x < 2Gj+i, the solution is given by 

tp'(x)\ _ f vfc \ /cos {k[x — Xj]) — sin (k[x — Xj]) 
i/)(x) J I -jr J ysin (k[x — Xj]) cos {k[x — Xj]) 




-) 

By recurrence, we obtain the solution for every x > in terms of a product of 
matrices. In particular, 

(3.2) ft? n+1 1)=^-i---A 1 f^ 1 -) 
where 
(3.3) 
and 




cos (kdj) — sin (k6j) 
sin (k9j) cos (k9j 




Thus, for instance, we see that a product of matrices of the form (1.6) corresponds 
to a generalised Kronig-Penney model of unit energy in which the o~Bj are double 
impurities. It is worth emphasising this point: between impurities, the Schrodinger 
operator itself produces the "rotation part" of the matrices in the product. There- 
fore, in order to associate a quantum model to the most general product of matrices, 
it is sufficient to use a potential made up of (suitably spaced) double impurities. 
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The case of negative energy, i.e. E = —k 2 , k > 0, is also of mathematical 



interest. Then Equation (3.2) holds with 




( . _ (Vk \ f cosh. (k6j) $mh(k6j) 

[ ' Aj -{0 ^) v sinh(fc^) cosh(fc^) 

3.2. The generalised Frisch Lloyd equation. It is physically reasonable to as- 
sume that the scatterers are randomly, independently and uniformly distributed. 
We denote by p the mean density of impurities. If we label the scatterers in order 
of increasing position along the positive semi-axis, so that Xj denotes the position 
of the jth impurity, then 

< X\ < X2 < x 3 < ■ ■ ■ 
and the spacings between consecutive scatterers are independent and have the same 
exponential distribution, i.e. 

(3.5) 6j ~ Exp(p) , p > . 
For this distribution of the 6j , 

n(x) := # {xj : Xj < x} 

is the familiar Poisson process. 

We shall be interested in the statistical behaviour of the Riccati variable 

(3.6) z{x) = . 
Its evolution is governed by 

(3.7) z' = -{z 2 + E) , xi { Xj }, 
and 

(3.8) z( Xj +) = (z(xj-)) , j G N . 

The "lack of memory" property of the exponential distribution (3.5 1 implies that 
the process {z(x)} thus defined is Markov. 



It should be clear from [1.1 and the previous subsection that, if we set k = 1, 



then the invariant measure v associated with the product (3.2) is precisely the 
stationary distribution of the Riccati variable. So we shall look for particular cases 
where this stationary distribution may be obtained in analytical form. 

To simplify matters, we also suppose in the first instance that the Bj are all the 
same, deterministic, and we drop the subscript. 

Let f(z;x) be the density of the distribution of the Riccati variable. Let h > 
and let dz denote an interval of infinitesimal length dz centered on the number z. 
Then 

f(z; x + h) dz = P (z(x + h) £ dz) 

CO 

= p ( z ( x + h) £dz n(x + h)- n(x) = i\ P (n(x + h) - n(x) = I) . 

1=0 

It is well-known (see [TB]) that, with an error of order o(h) as h — > 0+, 

fl-ph if£ = 
ph i(t=l 
if £ > 1 
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and so 



(3.9) f(z; x + h) dz = P (z(x + h) e dz n(x + h) - n(x) = OJ (1 - ph) 

+ P (^z(x + h) e dz n(x + h) - n(x) = 1 J ph + dz o(h) as h -t 0+ . 

The condition n(x + h) — n(x) = means that no Xj lies in (x,x + h), and so implies 
that the Riccati variable is governed solely by the differential equation (3.7) in this 
interval. Therefore, the first conditional probability on the right-hand side of (3.9) 
equals 

f(z+ [z 2 + E]h; x) [1 + 2zh + o(h)] dz as h -> 0+ . 

The condition n(x + h) — n(x) = 1 means that exactly one of the Xj lies in (x,x+h), 
and so the Riccati variable experiences a jump defined by Equation (3.8) in this 
interval. A simple calculation then yields, for the second conditional probability on 
the right-hand side of (3.9), 

P (z(x + h) edz n(x + h) - n(x) = lj = P {& (z(x ~)) G dz) + dz 0(h) 

= f (m- l {z);x) ^-^-{z)dz [I + 0(h)] as h -> 0+ . 
dz 



After reporting these results in Equation (3.9) and taking the limit as h —> 0+, we 
obtain a generalisation of the equation (6.69) in [31,, §6.7: 



d J x {z-x) = ^ z [{z 2 + E)f{z:x)]+p 



f(^- 1 (z);x) d ^(z)-f(z;x) 



The stationary distribution, denoted again / = f(z), therefore satisfies 
d 



[(z 2 + E)f(z)] +p 



0. 



More generally, if we permit the Bj to be independent random variables with a 
common distribution denoted k, then it is straightforward to derive the equation 



(3.10) -[(z 2 + E)f(z)]+p 



n(dB) 



SL(2,I 



f{®-\zj)^{z)-f{z) 



= 0. 



By integrating with respect to z, we obtain 



(3.11) 



(z 2 + E)f(z)+p 



K(dB) 



SL(2,] 



dt f(t) = N . 



The constant of integration N in this equation depends on E\ as will be ex- 
plained shortly, it represents the integrated density of states per unit length of 
the Schrodinger Hamiltonian for the potential (3.1) 19, 29, 31 . 



We shall refer to Equation (3.10), or to its integrated version (3.11), as the 
(generalised) Frisch-Lloyd equation. In the following sections, we shall consider 
again the particular point scatterers described in f|2j and exhibit choices of the 
measure K for which this equation can be converted to a differential equation. 
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3.3. The qualitative behaviour of the Riccati variable. It is instructive to 
think of the Riccati equation (3.7) in the absence of scatterers as an autonomous 
system describing the motion of a fictitious "particle" constrained to roll along the 
"potential" curve 

U(z)=Ez+^ 

in such a way that its "velocity" at "time" x and "position" z is given by the slope 
—U'(z); seeFigure[2] We may regard the occurence of the jumps in Equation (3.8) 
as a perturbation of this autonomous system, and the intensity p of the Poisson 
process as the perturbation parameter. 




z 



Figure 2. The "potential" U(z) associated with the unperturbed 
Riccati equation z' = —U'(z) = —(z 2 + E). 



Let us consider first the unperturbed system (i.e. p = 0). For E > 0, the 
system has no equilibrium point: the particle rolls down to — oo, and re-appears 
immediately at +oo, reflecting the fact that the solution i\> of the corresponding 
Schrodinger equation has a zero at the "time" x when the particle escapes to infinity. 
This behaviour of the Riccati variable indicates that every E > belongs to the 
spectrum of the Schrodinger operator. 



Equation (3.7) gives the "velocity" of the fictitious particle as a function of its 



position. Hence the "time" taken to go from 4-oo to — oo is 

dZ 7T 



z 2 + k 2 k 



On the other hand, the solution of the Frisch-Lloyd equation for E > and p = 
is the Cauchy density 

(3.12) /(z) = _JL_ with N= k _. 

Therefore the normalisation constant N may be interpreted as the reciprocal of the 
"time" that the particle takes to run through R. Another equivalent interpretation 
is as follows: recall that, when the particle escapes to — oo, it is immediately re- 
injected at +00 to commence a new journey through E. N may therefore also be 
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viewed as the current of the fictitious particle [31] , and the Rice formula 

lim z 2 f(z) = N 

z— >±oo 

can be understood as expressing a relation between the stationary distribution and 
a current of probability. This current equals the number of infinitudes of z(x) — 
i.e. the number of nodes of the wavefunction VK 2 -) — per unit length. By the 
familiar oscillation theorem of Sturm-Liouville theory, it is therefore the same as 
the integrated density of states per unit length of the corresponding Schrddinger 
Hamiltonian. 

By contrast, in the case E = — k 2 < 0, k > 0, the unperturbed system has an 
unstable equilibrium point at — k, and a stable equilibrium point at k. Unless the 
particle starts from a position on the left of the unstable equilibrium, it must tend 
asymptotically to the stable equilibrium point. The fact that the particle cannot 
reach infinity more than once indicates that the spectrum lies entirely in K + . The 
solution of the Frisch-Lloyd equation is 

f(z) = S(z - k) . 

Let us now consider how the occurence of jumps can affect the system. For E > 0, 



the jumps defined by Equation (3.8|, as long as they are finite and infrequent (i.e. 
p is small), cannot prevent the particle from visiting — oo repeatedly; the system 
should therefore behave in much the same way as in the unperturbed case, and 
we expect the density / to be Cauchy-like. In particular, the interpretation of the 
normalisation constant in terms of a probability current remains valid because, for 



z large enough, the deterministic part (3.7) of the evolution of the Riccati variable 



dominates the stochastic part (3.8). The situation for E < is more complicated. 
Roughly speaking, positive jumps, i.e. discontinuous increases of z, enable the 
particle to make excursions to the right of the stable equilibrium point z = k, 
but the particle can never overcome the infinite barrier and so it rolls back down 
towards k. On the other hand, negative jumps, i.e. discontinuous decreases of z, 
enable the particle to make excursions to the left of k. If the jump is large enough, 
the particle can overcome the potential barrier at — k and escape to — oo, raising 
the possibility that part of the spectrum of the Schrodinger operator lies in R_. 
For small p, we expect the density / to be large in the neighbourhood of z = k. 

We shall return to this useful particle analogy in later sections when we examine 
the detailed behaviour of the Riccati variable for specific random point scatterers. 



3.4. The reduced Lyapunov exponent. Knowing the density of the invariant 
measure, the calculation of the Lyapunov exponent reduces, thanks to Formula 
(1.3), to the evaluation of a multiple integral. We will show in this subsection that, 
if A is of the form (3.3) and the boundary matrix B is triangular, then this formula 
may be greatly simplified. As pointed out in S3.1 there is no loss of generality in 
setting E = k 2 = 1 since the parameter k may be re-introduced subsequently by 
rescaling. This has the advantage of making the calculation simpler. 

For definiteness, let us begin with the upper triangular case, i.e. &21 = 0. Then 
the density of the invariant measure satisfies the Frisch-Lloyd equation (3.10), again 
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with E = 1, and wc have 



.4 



Hence 



7m 



D 



(buz + b 12 ) 2 + (b 21 z + b 22 f 



(b 21 z + b 22 f [1 + 38(zf\ = b\ 2 [1 + @{zf 



Az 



n{dB) In 



A 



SL(2J 



2 lAz 



SL(2,1 



, x b\ 2 \l+M(z) 2 ] , x 

1 + 2^ 



and, after some re- arrangement, 

(3.13) 7/1 =E(ln|& 22 |) + i / dz / K(dB) ln[l + ^(z) 2 ] /(* 

* JR JSLf2,R) 



2 ' 1 



rc(cLB) ln(l + z 2 )/(z) 



SL(2,1 



Consider the second term on the right-hand side of the last equality: by changing 
the order of integration, and making the substitution y = 38 (z) in the inner integral, 
we obtain 



dz 



SL(2,: 



k(cLB) In [1 + S8{zf] f(z) 



\ f K(dB) f dzln[l + ^(z) 2 ]/(z) 

L JSL(2,R) JR 



SL(2J 



rullJ) j dylnfl + y 2 ]/^- 1 ^))^^) 



Next, we use the letter z instead of y, and change the order of integration again: 
Equation (3.13) becomes 



7/1 =E(ln|& 22 |) 
1 

+ 2 



dz 



n(dB) 



SL(2,1 



ln(l + z 2 ) 



Finally, by making use of Equation (3.10), and then integrating by parts, we arrive 
at the following formula: 



(3.14) 
where 
(3.15) 



7 ^ 7 



/OO 
dzzf(z), h 
-OO 



i=0. 
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This formula remains unchanged after restoring k by rescaling. A similar calculation 
may be carried out if, instead, B is lower triangular. Equation (3.14| then holds 
with 

(3.16) 7 :=j»E(ln|6ii|)-4 dz -f(z) , b 12 = . 

J -oo z 

The integrals in these expressions are Cauchy principal value integrals. 

We shall henceforth refer to 7 as the reduced Lyapunov exponent. Although our 
derivation of the relation between 7^ and 7 assumed that E > 0, we conjecture, 
on the basis of the numerical evidence obtained in all the examples we considered, 
that it holds also when E < 0. 

Such simplified formulae for the Lyapunov exponent are well-known in the physics 
literature [31] . The reduced Lyapunov exponent is the rate of growth of the solution 
of the Schrodinger equation: 



7 = lim — In \l ip{x) 2 + [ip'(x) 



IS! \ t '( .:■ I - -j- I i ' I X j 1 " 

x— >oo X 

Alternatively, using the stationarity of the process {z(x)}, 



and 



7= lim — In I "0 (ar) | + lim — In ^/l + z 2 (x) = lim — ln|^(x)| 

x— ¥00 x x—^oo x x^oo x 

7= lim — In + lim — In v/l/ ' z 2 (x) + 1 = lim — Inl^'fir)! 

x— foo X x^roo X x^roo X 



7 also provides a reasonable definition of (the reciprocal of) the localisation length 
of the system. 

The presence of the expectation term on the right-hand side of Equation (3.15) 
may, at first sight, surprise readers familiar with the case of delta scatterers, but 
its occurence in our more general context is easily explained as follows: between 
consecutive scatterers, the wave function is continuous and so, for x n < x < x n+1 , 
we can write 

(3.17) ln|^z)|=ln|V>(z„+)| + ^ dj/^lnfo&(|/)| 

= ln\ip(x n +)\ + f dyz(y). 

Let us denote by b^J the entry of B n in the ith row and jth column. If B is upper 
triangular, we have, at x n , 

1p(x n +) = &22 i>( x n-) 

and so the wavefunction is discontinuous there unless — 1- Reporting this in 
(3.17) and iterating, we obtain 

n(x) x 
ln|^(x)|=ln|V(0)| + Vln \b$ + / dy z(y) 

where n(x), as defined in S3.2 is the number of point scatterers in the interval 
[0, a;]. Using 

E(n(x)) = px 
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and the crgodicity of the Riccati variable, we recover by this other route the equa- 
tion (3.151 obtained earlier. It is now clear that the expectation term arises from 
the possible discontinuities of the wave function at the scatterers. To give two ex- 
amples: for the delta scatterer, the wave function is continuous everywhere, 622 = 1 
and so Equation (3.15 1 is just the familiar formula in [31 . For the supersymmetric 
scatterer, however, the wavefunction has discontinuitcs, 

622 = e~ w 

and so, as noted in [5T], the formula for the reduced Lyapunov exponent must 
include the additional term 

E(ln|6 22 |) - -E( W ). 
If B is lower triangular instead, it is more natural to work with 1/}' ': for x n < x < 



we have 



m|V'(*)| =m|V'(*n+ 
= m|V'(*n+)| 



dy — ln|V>'(y)| 



dy = ln|V>'Or„+)|-fc 2 



ln\^'(x n +)\-k 2 



dy 



dy 



z (y) 



Using the lower triangularity of B ni we obtain, at x n , 



By repeating our earlier argument, we recover Equation (3.16). 

3.5. Halperin's trick and the energy parameter. For the particular case of 
delta scatterers, Halperin [35] devised an ingenious method that, at least in some 
cases, by-passes the need for quadrature and yields analytical expressions for the 
reduced Lyapunov exponent. Let us give a brief outline of Halperin's trick and 
discuss some of its consequences. 

Halperin works with the Fourier transform of the invariant density: 



(3.18) 



F(x) 



/We- 



'd*. 



For the delta scatterer, the Frisch-Lloyd equation (3.11 1 in Fourier space is then 

E (e- ixu ) - 1 



(3.19) 



F"(x)-EF(x) -v 



F(x) 



-2ttNS(x). 



Let e > and integrate over the interval (—e,e). Using the fact that 



F'(-e)=-F'(s) 



and letting e — > 0+, we obtain 



N = Re[F'(0+)] 



Furthermore, since in this case b\\ = 1 and 621 = 0, Equation (3.151 leads to 

J — ( 



7 = 



:f(z) dz = -Im [F'(0+)] 
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These two formulae may be combined neatly by introducing the so-called char- 
acteristic function £1 associated with the system |32l 134] : 

(3.20) fl{E) := j(E) - mN(E) . 

Then Halperin observes that 

(3.21) 

where, with a slight abuse of notation, F is now the particular solution of the 



homogeneous version of Equation (3.191 satisfying the condition 



lim F(x) = . 

x— >+oo 

Thus the problem of evaluating 7 and N has been reduced to that of finding the 
recessive solution of a linear homogeneous differential equation. 



Equation (3.21) expresses a relationship between the density of states and the 
Lyapunov exponent — a relationship made more explicit in the Herbert-.] 'ones- 
Thouless formula [261 142] well-known in the theory of quantum disordered systems. 
A further consequence of the same equation is that, if the recessive solution F 
depends analytically on the energy parameter, so does the characteristic function. 
n should thus have an analytic continuation everywhere in the complex plane, save 
on the cut where the essential spectrum of the Schrodinger Hamiltonian lies. 

More generally, for an arbitrary scatterer, the Fourier transform F of the invari- 
ant density satisfies the equation 

(3.22) F"(x) - E F(x) - p n(dB) / dz f{z) = 0. 

isL(2,K) JR 1X 

We shall not make explicit use of this equation in what follows: instead, we shall 
obtain closed formulae for the characteristic function by making use of analytic con- 
tinuation. This is one important benefit of having retained the energy parameter. 

3.6. Random continued fractions. There is a close correspondence between 
products of 2 x 2 matrices and continued fractions. Let z be an arbitrary starting 



value, recall the definition (1.4) and set 

(3.23) z„, := Mi-i • • • M)(zo) ■ 

Then the sequence {z n } ne ^ is a Markov chain on the projective line, and v is fi- 
invariant if and only if it is a stationary distribution of this Markov chain. Now, 
reverse the order of the matrices in the product, set £0 = ^0 and 

(3.24) Cn :=^o°---°^n-i(Co)- 

Although, for every n, z n and £ ra have the same distribution, the large-n behaviour 
of a typical realisation of the sequence {z„}„ s n differs greatly from that of a typical 
realisation of the sequence {CnjneN [20] . 

{CrajnGN converges to a (random) limit, say C- Write 

. _ fa n b n 
n '~ \ c d 

Then 
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and so 

1/ci 



(3.26) C := ao/co , .... , 

dx/ci + a 2 /c 2 - —. , 

d 2 /c2 H 

It is immediately clear that, if A is /i-distributed, then 

1/c 2 
d/c + C 

has the same distribution as £. Hence the distribution of Q is /x-invariant. Further- 
more, if £ is independent of zq, then there can be only one /^-invariant measure. So 
/ is also the density of the infinite random continued fraction £. 

By contrast, behaves ergodically. Therefore the density / of the invari- 

ant measure v should be well approximated by a histogram of the z n . We have 
used this to verify the correctness of our results. 

4. Some explicit invariant measures 

4.1. Delta scatterers. In this section, we obtain invariant measures for products 
where the matrices are of the form 

_ (Vk \ (co S (ke) sin(fc<?)\ \ /I „ 



(4.2) A = 



'Vk \ /cosh(fc6») sinh(fc0)\ / ^ \ (\ u \ 
i J ^sinh(A:0) cosh(fc0) j Vifc J V° 1 ) 



where 9 ~ Exp(p) and u is a random variable, independent of 9, whose density 
we denote by g : M — >• M + . These products are associated with the generalised 
Kronig-Penney model for E = k 2 > and E = —k 2 respectively, in the case where 
(see Example [T]) 

*-c.i: 

The corresponding Frisch-Lloyd equation (3.11) is 

p rZ — U 

(z 2 + E)f(z)+p du dtf(t)g(u) = N. 

We change the order of integration; the equation becomes 

(4.3) N = (z 2 + E) f(z) + p [ dt K(z - t) /(f) 

JR 

where 

(4.4) K{x) 



— J g(u) du if x > 
e(u) du if x < 



Suppose that 
(4.5) ± u ~ Exp(g) . 
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We shall show that, in this case, the Frisch-Lloyd equation reduces to a first-order 
differential equation. For the sake of clarity, consider first the case u ~ Exp(g). For 
this choice of distribution, 



-e~i x if x > 
if x < 



So we have 



and equation (4.3) is 
N 



Differentiate this equation with respect to z: 
= ^ [(z 2 + E) f(z)] + pK(0+)f(z) + p 



K(x) 

K'(x) = -qK(x), K(0+) = -l, 
(z 2 + E)f(z)+p f dtK(z-t)f(t). 

dt K'(z-t)f(t) 



[(z 2 + E) f(z)} - pf(z) -qp dt K(z - t) f(t) 



= — [(z 2 + E) f(z)] - pf(z) -q[N-(z 2 + E) f(z)] 



This is the required differential equation. The case 
so we find, for the general case (4.51, 

d 



Exp(q) is analogous, and 



(4.6) 



dz 



[(z 2 + E) f(z)] - p f(z) ± q [(z 2 + E) f(z)] = ±qN . 



We seek the particular solution that satisfies the normalisation condition 



(4.7) 



f(z)dz=l 



This condition fixes the constant of integration N, and hence provides an expression 
for the integrated density of states for the Schrodinger Hamiltonian. 



4.1.1. Product of the form (4-1). For E = k 2 > 0, this leads to 
±qN 



(4.8) f(z) 



z 2 + k 2 



exp 



Tqz 



p z 
— arctan — 
k k 



exp 



=poo 



p t 

±qt arctan — 

k k 



dt. 



The density of the Riccati variable is plotted in Figure [3] for (a) positive u n 
and (b) negative u n . The continuous black curves correspond to a low density of 
scatterers (small p, compared to k and 1/q) and are reminiscent of the Cauchy law 
obtained in the absence of scatterers. Recall that the effect of the nth scatterer on 
the Riccati variable is described by the equation 

z(x n + ) = Z(x n -) + U n ■ 

When the u n are positive, any increase in the concentration of the scatterers pro- 
duces a decrease in the current and so the distribution is pushed to the right; see 
the blue dashed curve in Figure [3] (a). On the other hand, when the u n are nega- 
tive, any increase in the concentration of the scatterers leads to an increase in the 
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current of the Riccati variable and so spreads the distribution; see the blue dashed 
curve in Figure [3] (b). 
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(b) 



Figure 3. Plots of the invariant density / for delta scatterers and 
positive energy E — k 2 = +1. Black continuous lines correspond 
to a low density of scatterers, i.e. ~ Exp(p) with p = 1/4, and 
blue dashed lines to a high density, i.e. p = 4. (a) u ~ Exp(l); (b) 
— u ~ Exp(l). 



4.1.2. Product of the form (4-2). For E = —k 2 < and u ~ Exp(g), one must take 
iV = to obtain a normalisable solution — a reflection of the fact that the essential 
spectrum of the Schrodinger Hamiltonian is M + . Then 



(4.9) 



/(*) = c 



k 2 



- k 
~+k 



l(fe,oo)(^) 



where C is the normalisation constant. 

By contrast, in the case —u ~ Exp(g), one finds 



(4.10) 



where 



/(*) 



9 7V 



-k 2 



z — k 


P 

2k 


r c(z) 

/ e-«* 

J z 


t- 


-k 


p 

2 k 


z + k 




t - 


- k 





dt 



c(z) 



DO 



if z > fc 
if z < k 

The invariant density / is plotted in Figure [I] for the cases (a) u ~ Exp(g) and 
(b) — u ~ Exp(g) respectively. The shape of the distribution can again be explained 
by using the qualitative picture of §3.3| For positive u n , the sharp peak obtained 
for a small concentration p of scatterers (black continuous line) reflects the trapping 
of the process {z(x)} by the potential well at z — k; recall Figure [2j When the 
concentration of scatterers is increased, the Riccati variable experiences positive 
jumps more frequently and so the distribution spreads to the right (blue dashed 
curve). For negative u n , the jumps can take arbitrary negative values. This enables 
the "particle" to overcome the barrier at z = —k, and so we have a non-zero current 
N (i.e. a non-zero density of states). This effect is enhanced as the density of the 
scatterers is increased (blue dashed curve). 
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(a) (b) 

Figure 4. Plots of the invariant density / for delta scatterers and 
negative energy E = —k 2 = — 1. Black continuous lines correspond 
to a "low" density of scatterers, i.e. 9 ~ Exp(p) with p = 1/4, and 
blue dashed lines to a "high" density p = 4. (a) u ~ Exp(l); (b) 
— it Exp(l). 



4.1.3. Calculation of the characteristic function. Using the invariant measure, it is 
trivial to express the integrated density of states and the Lyapunov exponent in 
integral form. Such integral expressions are particularly useful when studying the 
asymptotics of N and 7 in various limits. Even so, it is worth seeking analytical 
expressions (in terms of special functions) for these quantities, as they sometimes 
reveal unexpected connections to other problems. 

Recalling the discussion in |3.5| we begin with a straightforward application of 
Halperin's trick. For ±ti ~ Exp(g), we have 

E (e- iux " ' 



1 ± ix/q 

and so the homogeneous version of Equation ( 3.19[ ) is 

P 



F"(x) + 
The recessive solution is 

F(x) = 



-E 



±q 



F(x) = 0. 



W 



{2y/^E [ix ± q}) 



W 



{±2V=Eq) 



where W a ,b is the Whittaker function Q}[22|. We deduce that, for E outside the 
essential spectrum of the Schrodinger Hamiltonian, 

W ^ , (±2V^Eq) 

(4.11) Q(E):=j(E)-mN(E) = -2V^E: 



.4 (±2V=^g) 



This formula for the characteristic function was discovered by Nieuwenhuizen 
In particular, for k real, 

W'^ , (T2ikq) 

7 (*») - far^) = n(k 2 + 10+) = 



22 



ALAIN COMTET, CHRISTOPHE TEXIER, AND YVES TOURIGNY 



and 



7 (-/s 2 ) - mN(-k 2 ) = n{-k 2 + i0+) = -2k- 



WLt i (±2fcg) 



W- P i (±2fcg) ' 

2fc I 2 

In the case u ~ Exp(g) there is an alternative derivation of this formula which 
does not require the solution of a differential equation: start with the explicit form 
of the invariant density / for E — —k 2 < 0, given by Equation (4.9 1. By using 
Formula 3 in |23| . §3.384, we obtain the following expression for the normalisation 
constant: 



C := 



z 2 — k 2 



- k 
Tk 



2k 



(I) W *,\ « 



The reduced Lyapunov exponent 7 may then be obtained easily by noticing that 
differentiation with respect to the parameter q yields an additional factor of z in 
the integrand. Hence, for k real, we find 

,00 q W' =s _ i (2kq) 

7 (-fc 2 ) = / dzzf(z) = -— InC* = -2k 



dq 



W-p 1 (2kq) 



Since N = for E < 0, this yields 
(4.12) 



W'^ , (2kq) 

n(-k 2 ) = -2k — — - 

V ' W-v 1 (2kq) 

2fc ' 2 



Now, the half-line E < lies outside the essential spectrum of the Schrodinger 
Hamiltonian because N = along it. Hence fl is analytic along this half-line, 



and we see that the "+ c ase" o f our earlier Equation (4.11) is simply the analytic 
continuation of Equation (4.12). In particular, the formula in the case E — k 2 > 
may be deduced from the formula in the case E = — k 2 < by applying the simple 
substitution 

k 1— > —ik . 

4.2. Delta— prime scatterers. Products of matrices of the form 

Vfc \ A-osf/.-fl] -,m(kt)\\ ( -V \ / I 1 ) N 



(4.13) 



.4 



(4.14) 



A = 




cos(fc#) — sin(fc0) 
sin(fe#) cos(fc#) 




cosh(fc<9) sinh(fc0)\ 
sinh(fc6>) cosh(fc0)i 1 \/k 




where 9 ~ Exp(p) and v is a random variable independent of 9, are associated with 
the delta-prime scatterer (see Example [2| 

'1 0^ 

.v 1, 



B = 



The Frisch-Lloyd equation (3.11) for this scatterer is 
(4.15) (z 2 + E)f(z)+p [ dv j 1 '" dtf(t)g(v) = N 



where g is the density of v. The calculation of the invariant measure in this case 
can be reduced to the calculation of the invariant measure for some Kronig-Penney 
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model with delta scatterers. For instance, in the negative energy case (4.141 with 
k = 1, we have 



.4 



1 

1 



/ cosh 9 sinh 9 
lsinh6> cosh# 



The similarity transformation of the matrix A on the left corresponds to the trans- 
formation z i — y X/z of the Riccati variable. So the invariant densities for the delta 
and the delta-prime cases are in a reciprocal relationship. Accordingly, replace z 



by l/z in Equation (4.15) and set 



1 



g(z) = -zf(l/z) 



Then 

(4.16) N = (1 + Ez 2 ) g(z) - p I dv 



dtg(t)g{v) 
(l + Ez 2 )g(z)-p 



dtK(z - t) g(t) 



where K is the kernel defined by Equation (4.4 1. This equation for g is essentially 



the same as Equation (4.3 ) save for the sign of p and the dependence on the energy. 
For the case 

±v ~ Exp(g) 

this equation can, by using the same tricks as before, be converted into a differential 
equation which is easy to solve. 



4.2.1. Product of the form (4-13). For E = k 2 > 0, the upshot is 



(4-17) f{z) 



±qN 
z 2 + k 2 



exp 



p k 
— arctan — 

k z 



l/z 



exp 



=poo 



±qt H — arctan(fci) 
k 



dt. 



Plots of the distribution are shown in Figure[5]for (a) positive v n and (b) negative 
v n . These plots differ somewhat from those obtained in the case of delta scatterers, 



and we can use the particle analogy of S3.3 to explain the differences. The jump of 
the particle associated with the nth delta-prime scatterer is given implicitly by 

(4.i8) -r L ^ = ^—^ +Vn - 

z(x n +) z(x n -) 

The strongly asymmetric distribution obtained for negative v n (Part (b) of Figure 
[5]) can be explained as follows: starting from +00, the particle experiences its first 
jump at "time" x%, and its value after the jump is approximately 1/vi < 0. In fact, 
for z negative and small in modulus, the invariant density resembles very closely 
that of 1/vi, i.e. 



/(*) 



as z — > 0- 



Thereafter, the particle proceeds towards —00. In particular, if p is large, then 
the expected value of X\ is small, and the particle spends hardly any time on the 
positive semi-axis. 
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(a) 



(b) 



Figure 5. Plots of the invariant density / for delta-prime scat- 
terers and positive energy E = k 2 = +1. Black continuous lines 
correspond to a "low" density of scatterers, i.e. 9 ~ Exp(p) with 
p = 1/4, and blue dashed lines to a "high" density p = 4. (a) 
v ~ Exp(l). (b) -v ~ Exp(l). 



Product of the form (^.l^): for E = —k 2 < and v ~ Exp(q), we find 

(4.19) /<=)^-'^(^) "«« 
where 

(4.20) c = i T ii) w -^J 2 ^- 

When E = —k 2 < and — v ~ Exp(g), we obtain 



(4.21) 
where 



z 2 _ fc2 



c(z) 



2 — k 


V 
2k 


r l/z 

/ e-«* 


1 + kt 


z + k 




7 C ( 2 ) 


l-kt 



dt 



00 if < z < l/k 
— k otherwise 



Again, we can try to understand the qualitative features of the density function 



/ for E < by invoking the particle analogy of [3.3 In view of Equation (4.18), 
when v n ~ Exp(g) and z(x n — ) > 0, the value of the Riccati variable decreases 
but can never become negative. So the particle, once it passes to the left of the 
equilibrium point at z = k, must remain trapped there. This explains why the 
density is supported on (0, fc); see Figure [6] (a). By contrast, when — v n ~ Exp(g), 
the jumps are unrestricted; the "particle" can escape over the potential barrier at 
z = —k infinitely often, leading to a non-zero current and a density / spread over 
K. This is shown in Figure [6] (b) . 



4.2.2. Calculation of the characteristic function. We begin with the case v 
and E = -k 2 < 0. 



Exp(g) 



The invariant density is then given by Equation (4.19). Using Formula (3.16) 



for the reduced Lyapunov exponent and the expression (4.20) for the normalisation 
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-0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 ' -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 



(a) (b) 

Figure 6. Plots of the invariant density / for delta-prime scat- 
terers and negative energy E = — k 2 = — 1. Black continuous lines 
correspond to a "low" density of scatterers, i.e. 9 ~ Exp(p) with 
p = 1/4, and blue dashed lines to a "high" density p = 4. (a) 
v ~ Exp(l). (b) -v ~ Exp(l). 



constant C, we find 

g W^, (2q/k) 

,y ' dq W^p i (2q k) 

2k > 2 

Since N(—k 2 ) = in this case, analytic continuation yields 



(4.22) fi(J5) = -2V-E 



The characteristic function in the case — u ~ Exp(g) is the same, except that g 
becomes — g. In particular, for E — k 2 > and ±i> ~ Exp(g), we obtain 

WL± , (±2iq/k) 
7 (*») - i.A^ 2 ) = ^ 2 4- 10+) = 21 ^-- • 

2fc ' 2 

An alternative derivation of these results could use the correspondence between 
the delta and delta-prime cases alluded to earlier. 



4.3. Supersymmetric scatterers. We now consider products where the matrices 
are of the form 

_ (Vk \ /cos(fcfl) - sin(fc0)\ 0\(e» \ 

1 ' A ~\0 \sm{k6) cos(fcfl) )[ Vk) \ e~™ J 

or 

_ /V* \ fcosh(fcfl) sinh(fcfl)\ M 

1 ' A -[0 ^Msinh(fe0) cosh(*0)Ho Vk \0 e~ w J 
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where 9 ~ Exp(p) and w is a random variable independent of 6. These products 
arise in the solution of the generalised Kronig-Penney model with the supersym- 
metric interaction of Example [3j i.e. 





B = 







Let g denote the density of w. The Frisch-Lloyd equation (3.11 1 is 
(z 2 + E)f(z)+p [ dw [ dtf(t)g(w)=N. 



After changing the order of integration, this becomes 

f°° (\ z 

(4.25) N = (z 2 + E)f(z)+p K (-In- 



fit) dt 



where K is the kernel defined by Equation (4.4 1. 
Let 

±u> ~ Exp(g) . 

Then the kernel is supported on R± and satisfies the differential equation 
K'{x) = TqK(x) , K(0±)=Tl- 

We deduce 

(4-26) A [ {z 2 + E) /(z) ] _ p f(z) ± q z l+A f{z) = ±y z n 

where N is the integrated density of states. 

Before going on to solve this equation, let us make a general remark: in the 
supersymmetric case, if one knows the invariant density, say /+ , for a certain dis- 
tribution of the strength w, then one can easily deduce the invariant density, say 
/_, when the sign of the strength is reversed. For instance, in the case E = 1, the 
relationship between /_ and /+ is simply 



This relationship can be deduced directly from the form of the matrices in the 
product (4.23). It is also connected with the fact that changing the sign of the 
superpotential W in Example [3] corresponds to swapping the functions <j> and ip — 
a manifestation of the so-called supersymmetry of the Hamiltonian. 



4.3.1. Product of the form (4-23). For E = k > and ±w ~ Exp(g), we have 



(4.27) f(z)=N- 



where 



z 2 + k 2 



exp 



p z 
— arctan — 
k k 



c±(z) 



c+(z) = and C—(z) 



i 2 exp 



oo if z > 
— oo if z < 



P arctan * 
k k 



dt 

T 



For the supersymmetric scatterer, 
(4.28) m(z) 



e 2w z. 
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Hence, for w ~ Exp(q), the jumps increase the Riccati variable if it is already pos- 
itive, and decrease it otherwise. Furthermore there is no bound on the magnitude 
of the jumps. It follows that the effect of increasing the density p of the scatterers 
is to decrease the density / on K_ , and to increase it on R + . This is in agreement 
with the plots shown in Figure [7] (a) . 

For — w ~ Exp(q), we observe the opposite effect: as shown in Figure [7] (b), for 
increasing p, the density / is lowered on R + and raised on K_. The asymmetry 
of the plots for a negative strength w is readily explained by using the particle 
analogy: starting at +00, the particle rolls down the potential, spurred along by 
the impurities, and quickly reaches the origin. Once the particle crosses over to the 
left, the impurities work against the downward force and tend to push the particle 
back towards the origin. 




z 

(a) 



z 

(b) 



Figure 7. Plots of the invariant density / for supersymmetric 
scatterers and positive energy E = k 2 = +1. Black continuous 
lines correspond to a "low" density of scatterers, i.e. 6 ~ Exp(p) 
with p = 1/2, and blue dashed lines to a "high" density p = 2. (a) 
w ~ Exp(l). (b) — w ~ Exp(l). 



4.3.2. Product of the form (4-24)- For E = —k 2 < and w ~ Exp(g), we must 
take N = in Equation (4.26) to obtain a normalisable solution. This is consis- 
tent with the well-known fact that the spectrum of a supersymmetric Schrodinger 
Hamiltonian must be contained in R + . Hence 



(4.29) 



z 2 — k 2 

For — w ~ Exp(q), the solution is, instead, 



(4.30) 

By Formula 8 in 



/(*) = CI 1 



Z2 



k 2 



z + k 

k — z 
k + z 



§3.197, 



(4.31) 



C+ = fc T 5 



2-1 



2k 



P_ 
2k 



l(fc,oo)(^) • 



1 (o,k)(z) ■ 



1, 



2k 



l;-l 



where B is the beta function and 2-F1 is Gauss's hypergeometric function. 
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(a) 



(b) 



Figure 8. Plots of the invariant density / for supersymmctric 
scatterers and negative energy E = — k 2 = — 1. Black continuous 
lines correspond to a "low" density of scatterers, i.e. 9 ~ Exp(p) 
with p = 1/4, and blue dashed lines to a "high" density p — 4. (a) 
w ~ Exp(l). (b) —w ~ Exp(l). 



Plots of the invariant density are shown in Figure [8j As before, the particle 
analogy helps to explain their qualitative features: in view of Equation (4.28 ), when 
w ~ Exp(q), the "particle" must eventually end up to the right of the equilibrium 
point z — k; see Figure [8] (b). By contrast, when — w G Exp(g) and z > 0, the 
Riccati variable remains positive but its value decreases at every jump. Hence, in 
this case, the support of the invariant density is (0, fc); see Figure [8] (b). 



4.3.3. Calculation of the characteristic function. The essential spectrum — and 
hence also the characteristic function — is invariant under a change of sign of the 
strength w. So we need only consider the case w ~ Exp(g). For E = —k 2 , we find, 



by using Equation (4.291 



/oo 
dzzf(z) 
-OO 



CI 



-l+l 



dz 



k 2 



The normalisation constant C+ is given explicitly by Formula (4.31 1, and a similar 
formula is available for the definite integral; it suffices to replace q/2 by q/2 — 1. 
Using (3.15) and the fact that N{—k 2 ) = 0, the result is 



Q(-k 2 



P -+k 



2 k 



1 i- JL 

2' 2fe 



2. 

2 ' 



-1 



B (P_ 1 
D \ 2k ' 2 



i)^i(& + i>f + i 



2fc 



1;-1 



This formula extends to other values of the energy by analytic continuation; it 
suffices to replace k by \J—E. In particular, for E = k 2 > 0, the characteristic 
function is obtained by replacing k by — ifc. The density of states and the Lyapunov 



exponent may then be deduced from the formulae (see Equation (3.20)) 



N 



1 



Im Q, and 7 = Re 51 . 
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5. Extensions 

In this final section, we consider possible extensions of our results: (1) to another 
scatterer; (2) to another distribution of the strength of the scatterers and (3) to 
another distribution of the spacing between consecutive scatterers. 



5.1. Double impurities. The decomposition formula (1.6 1 gives a formal cor- 
respondence between products of 2 x 2 matrices and generalised Kronig-Penney 
models of unit energy where the point scatterers are double impurities. In the 
particular case 

6 ~ Exp(p) , 

the density of the invariant measure solves the Frisch-Lloyd equation 



(5.1) N = (z z + l)f(z)+p dudw e (u,w) f(y)dy. 

Jrjr J z 

We have already considered the cases where to vanishes almost surely (the delta 
scatterer) or u vanishes almost surely (the supersymmetric scatterer) . The purpose 
of this subsection is to consider the truly multivariate case where u and w are 
independent and 

u ~ Exp(<j d ) , w ~ Exp(q s ) . 

We shall show that the corresponding invariant density / solves the differential 
equation 

(5.2) — [2z (if' ~pf)} + (q s + 2q d z) {if 1 -pf) + q d q s p = q d q s N 
az 

where 

V := (z 2 + l)f(z) 



and N is independent of z. To derive this equation from (5.1), we shall consider 
the cases z < and z > separately. 
Consider the latter case; we write 

g{u,w) = Q d {u) 0g (to) . 

By changing the order of integration, we find 



dto0 s (to) / f{y)dy 

= - f f(y) dy + f U K. (\ In f(y) dy , 

Jz-u J-u V 2 y + uj 



where 

' - e - q ° x ifx>0 



K s {x) 



if x < 



So Equation (5.1) becomes 



N = cp(z)-p dug d (u) f(y)dy 



Z — U 



dug d (u) I K s ( - In ) f(y) dy . 

2 y + u 
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and, by using integration by parts for the first integral on the right-hand side, we 
hnd 



(5.3) N = <p(z)-p e- qdU f(z-u)du 



+ P dug d (u) 



h\ I Jin — 

.2 y + 



f(y) dy- 



Next, differentiate this equation with respect to z. By exploiting the identities 

K q s K s , K s (0+) = 1 

and 



du 



~ qdU = -q d e~ qdU . 



we deduce 

/•OO 

(5.4) 2z W(z) - V f(z)] + q s [<p(z) — N] = q sP / e"^" f(z - u) du . 

Jo 

The integral term may be eliminated by differentiating once more with respect to 
z, and we obtain eventually Equation (5.2). 

The same equation is obtained if, instead, z < 0. It is a trivial exercise to adapt 
these arguments to cater for cases where one or both of u and w is always negative. 
We do not know how to express the solution of this second-order linear differential 
equation in terms of known functions, except in the limiting cases 



q d fixed , q s )■ oo , 



and 



q d -> oo , q s fixed , 
that have already been considered in |4.1| and |4.3| respectively. 



5.2. Delta scatterers with a gamma distribution. The equation (4.3) with the 
kernel (4.4) can be reduced to a purely differential form whenever g solves a linear 
differential equation with piecewise constant coefficients. For instance, suppose that 

±u ~ Gamma(2, 1/q) , 

i.e. 

q(u) = ±q 2 ue^ u l R± {u) . 

Then 

d 



±qj g = 0. 

Using the same trick as before, we obtain the following differential equation for 
cp:=(z 2 + E)f: 

2 



(5.5) 



±q\ cp-p 



dz 



±2q 



z 2 + E 



= q 2 N . 



Suppose that E — k 2 > and use the ansatz 



ip(z) := exp 



p z 
^fqz + — arctan — 
k k 



h(z). 



Then 
(5.6) 



(z 2 + k 2 ) h" + pti T pqh = . 
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This equation may be solved in terms of hypergeometric series and, by imposing 
suitable auxiliary conditions, we can find a particular solution h that is positive. 
The method of variation of constants then yields 



(5.7) /(*) = 



where 



q 2 N 
j2 + p 



exp 



p z 
^fqz + — arctan — 
k k 



h(z) 
exp 



=poo 



p t 

±qt — — arctan — 
k k 



H(t) dt . 



H(z) 



h 2 (z) 



e ±qt h(t)dt. 



=poo 



We now return to the calculation of the function h{z) appearing in this formula. 
The general solution of Equation (5.6) takes a remarkably simple form when 



Indeed, substitute 
(5.8) 



K z ) : = T 2 ' ' a ° = 1 



i=0 



into Equation (5.6|. This yields a recurrence relation for the af. 

k 2 a l+2 = -pa i+1 + \pq-i(i- 1)] aj_i = , i = 0,l, ... . 

By choosing a\ so that dj+i = 0, the infinite series reduces to a polynomial, say 
Pi i 



P 2 (z) = 1 + 2 



pz 



+ 2 



(p 2 +4k 2 )z 
P 3 (z) = 1 + 6 -^-s ^Awv + 18 



p 2 + 2k 2 p 2 + 2k 2 

2 

24 



p(p 2 + lQk 2 ) p 2 + 10k 2 p(p 2 + 10k 2 ) 

pip 2 + 16 k 2 )z (p 2 + 6k 2 )z 2 

Ay ' p 4 + 28 P 2 k 2 + 72k 4 p 4 + 28 p 2 k 2 + 72k 4 



240 



pz 



360 



p 4 + 28p 2 fc 2 + 72fc 4 p 4 + 28p 2 /c 2 + 72fc 4 

Another solution may be found by setting 



h(z) := (z 2 + fc 2 )exp 



p z 
— — arctan — 
fe fc 



££a & o = 1 ' 



4=0 



Then 



A; 2 & i+2 =p6 l+1 + [pg- (i + 2)(i + !)]&< = 0, i = 0, 1, 



By choosing b\ so that bj-\ — 0, this series reduces to another polynomial, say Qj 

Q 2 (z) = 1. 



Q 3 {z) = l--z. 

P 

pz Z?" 
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Hence the general solution of Equation (5.6| is 
(5.9) h(z) = aP s (z) + c 2 Q 3 {z){z 2 + fc 2 )exp 



p z 
- — arctan — 
k k 



Even with such detailed knowledge, it is not straightforward to identify the 
particular solution h that yields the density. We end with the remark that the 
characteristic function may, nevertheless, be constructed by using Halperin's trick: 
in this case, the homogeneous version of Equation ( |3.19 ) is 

The solutions are expressible in terms of Whittakcr functions; in particular, for k 

n{k 2 ) = 1 (k 2 ) - mN{k 2 ) = 2ik — 2k ' 2 . — . . 

V ' n ' V 1 W. P (T2ifcg) 

1 2k > 2 

This result was originally found by Nieuwenhuizen |35j . 

5.3. An alternative derivation of the Frisch Lloyd equation. In deriving 
Equation (3.11), we made explicit use of the fact that, when the spacing 8j :— 
Xj + i — Xj is exponentially distributed, 

n(x) := # {xj : Xj < x} 

is a Poisson process. In this subsection, we outline an alternative derivation of the 
Frisch-Lloyd equation which generalises to other distributions of the 9j. 

There is no real loss of generality in setting E = 1. We use the decomposition 

( cos — sin 9\ ^ 
ysinfl cos 8 J 

Then 
where 

. . z cos 8 — sin 8 
z sm 8 + cos 8 

Denote by g the density of the random variable 8 and by k the distribution of B. 
Equation (1.5| for the invariant density / then becomes 

f{x) = f d9 g(9) f K (dB) [f o ^ o o <%_ e ] (z) 

JR JSL(2,R) dz 

f (\ c %~ 1 Bid 

d9g(8) K(dB)[fo^- 1 ](w) a ^—(w) C p 

JSL(2,R) az az 

where 

_ , . z cos 8 + sin 8 

w := M- e {z) = — — . 

—z sm 9 + cos 9 

The same equation can also be written in the more compact form 
(5.10) f(z)= / d8g(8)— I K(dB) / dtf(t). 



Now, 



SL(2,R) JO 



dw . o. dw 
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Hence, if we multliply Equation (5.10) by 1 + z , we obtain 



(l + z')f(z) = / d9g(0) 
Jn 

Next, differentiate with respect to z: 
<l [(l + z 2 )/(z)] = f d9g{9) °' 



86 



SL(2,1 



k(cLB) 



At fit) 



o 



dz 



dzdO 



K(dB) I dtf(t) 



SL(2J 



d9g{9) 



d 2 
d9dz 



SL(2,1 



K{dB) / dtf{t) 



We may then use integration by parts for the outer integral; in the particular case 

g(9)^ P e-P 9 l [0tOo) 



the Frisch-Lloyd equation (3.11) follows easily after invoking (5.10) once more. 



We can use the same trick whenever the density of 9 satisfies a linear differential 
equation with constant coefficients. For instance, in the case 

9 ~ Gamma(2,p) 

it may be shown that 



dz 
where 



[(1 + z VM] ~ WW + P 2 f(z) = P 2 / K(dfl) [/ ■ 

iSL(2,R) 



d^\ 



V {z) = (l + z 2 )f{z). 
6. Conclusion 

In this article we have studied the invariant measure of products of random ma- 
trices in SL (2, E). This study relied on the correspondence between such products 
and a certain class of random Schrodinger equations in which the potential consists 
of point scatterers. We have considered several instances of this correspondence: 
delta, delta-prime and supersymmetric scatterers. By generalising the approach 
developed by Frisch & Lloyd for delta scatterers, we have obtained an integral 
equation for the invariant density of a Riccati variable; this density yields the in- 
variant measure of the product of random matrices. For the three cases of point 
scatterers we have obtained explicit formulae for the invariant measures. These are 
the main new results of this paper. 

The integrated density of states and the Lyapunov exponent of these models 
were also calculated. Two approaches were used for this purpose: the first is 



"Halperin's trick" and is specific to the case of delta scatterers (cf. section 3.5); 
the second uses analytic continuation of the characteristic function and depends 
on the explicit knowledge of the invariant measure in some interval of the energy 
outside the spectrum. By the first of these methods we have recovered the results 
of Nieuwenhuizen in the case of delta scatterers. By the second method we have 
found new explicit formulae for the integrated density of states and for the Lyapunov 
exponent in the cases of delta-prime and of supersymmetric scatterers. 

All these analytical results were obtained when the spacing between consecutive 
scatterers, as well as the impurity strength, have exponential distributions. Possible 
extensions to the gamma distribution were also discussed. 
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A more complicated type of scatterer, combining the delta and the supersym- 
metric scatterers, has also been examined. We called this scatterer the "double 
impurity" ; it is interesting because every product of matrices in SL (2, M) may in 
principle be studied by considering a Schrodinger problem whose potential consists 
of double impurities. Although we succeeded in deriving a differential equation for 
the invariant measure associated with a particular distribution of such scatterers, 
we were unable to express its solution in terms of known functions. 

In this paper we have played down the physical aspects of the models. Apart from 
the inverse localisation length and the density of states, there are other physical 
quantities that bear some relation to the Riccati variable and whose statistical 
properties are of interest. Let us mention three of them: the most obvious is 
the phase of the reflexion coefficient on the disordered region; for a semi-infinite 
disordered region, its distribution is trivially related to the invariant density of 
the Riccati variable [SJ 122] ■ Another quantity is the Wigner time delay (the 
derivative of the phase shift with respect to the energy) ; it has been considered in the 
contexts of the Schrodinger [28, 40 and Dirac 39J equations. A third quantity is the 
transmission coefficient (i.e. conductance) [H 1 16} 15?]. The study of the distributions 
of the Wigner time delay and of the transmission coefficient is mathematically more 
challenging because it requires the analysis of some joint distributions; for this 
reason it has been confined so far to limiting cases. 

Some of the physical aspects arising from our results will form the basis of future 
work. 



Appendix A. Scattering, transfer and boundary matrices 

We discuss in this section the relationship between the scattering matrix S, 
defined by (2.9), and the boundary matrix B, defined by (2.3). Here we need not 
assume that the scatterer is necessarily pointlike; the scattering matrix and the 
corresponding boundary matrix could equally well describe the effect of a potential 
supported on an interval. 

We first write the scattering matrix in terms of transmission and reflexion prob- 
ability amplitudes t, t' and r, r' : 



S 



r t' 
t r' 



Current conservation implies 



| a out|2 + | a out|2 = | fl in|2 + | ffl in|2 

and so forces the scattering matrix to be unitary, i.e. S <E U(2). The constraints 
on the coefficients, namely 

|rf + \t\ 2 = \r'\ 2 + \t'\ 2 = |rf + \t\ 2 = \r\ 2 + \t'\ 2 = 1 , 

r'/t' = -r/i and r/t' = -VJt, 

are conveniently built into the following paramctrisation, which also illustrates the 
factorisation [7(2) = U(l) x SU(2): 

(A.1) S = ie»( e "^7 -'^V 

\ -VlfXy/T e-"fy/T^Tj 
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This representation of the scattering matrix is interesting because the four real 
parameters have a clear physical interpretation: r G [0, 1] is the probability of 
transmission through the scatterer; 9 is the global phase of the matrix, i.e. 

detS = -e 2if> . 

It is sometimes referred to as the "Friedel phase" since it is the phase appearing in 
the Krein-Friedel sum rule relating the local density of states of the scattering region 
to a scattering property. The phase ip is a measure of the left-right asymmetry 
(tp = or 7r corresponds to a scattering invariant under x — > —x). Finally the 
phase x is of magnetic origin, since time reversal corresponds to transposition of 
the scattering matrix. 

Next, we introduce the transfer matrix T relating left and right amplitudes: 

/ 1/i -r7t\ 

where T — 




-r/t' 

This matrix is useful when considering the cumulative effect of many scatterers 
because it follows a simple composition law. Again, current conservation implies 
that the transfer matrix is unitary: 

|a° ut | 2 - |a^| 2 = \a™\ 2 - |a° ut | 2 . 

In other words, T £ {7(1, 1) (note that detT = t/t' = e 2ix ). 

The boundary matrix is also a "transfer matrix" in the sense that it connects 
properties of the wavefunction on both sides of the scatterer. The relation between 
T and B is easily found: from 



W(°±)y v 1 




we deduce 

B = UTU^ where U 



ix/4 fik -ik 



2k VI 1 



of Equation ( 2.3 ): 



Then, using the parametrisation (A.ll, we arrive at the following alternative form 



c lx f cos 9 — sin ip \J\ — t —k [sin 9 + cos tp y/l — r] \ 

(A.2) B = —= . r . I 

v r V t L sm " ~ cos f v 1 — T \ cos 6* + sin t/3 V 1 ~ T / 

In particular, this expression shows clearly that 

e~ ix B G SL(2,R) . 

In one dimension, if a magnetic field is present, it may always be removed by a 
gauge transformation. Furthermore, setting the magnetic phase in the exponential 
factor to zero does not affect the spectrum of the Schrodinger operator. Hence 
there is no loss of generality in restricting our attention to the case B € SL(2,R). 

We end this appendix with some examples of scatterers, expressed in terms of 
the parameters x> T i @ an d <p- 

Example 1. For r = 1, ip = and \ — 0, B is the matrix describing a rotation 
of angle 9 — kt. In this case, the "scattering" is equivalent to free propagation 
through an interval of length £. 
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Example 2. The scattering matrix for the delta impurity may be written as 

10 /isin0 cos#" 



where 



cos 9 i sin 9 



= — arctan — e (— 7r/2, tt/2) . 



The other parameters are given by x = 0, 

fo ifu<0 



V = 



7T if u > 



and 



1 + 



2/r 



Example 3. For the delta-prime scatterer, the scattering matrix S has the same 
form as in the previous example, but this time with 

9 = arctan ^ e (-tt/2, tt/2) , 




and 



Example 4. The supersymmetric scatterer corresponds to taking \ = 6 = 0, 
99 = — 7r/2 and r = sech 2 u>. 
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